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Formulas are presented for the recursive generation of four-body integrals in which the integrand 
consists of arbitrary integer powers (> —1) of all the interparticle distances Tij, multiplied by an 
exponential containing an arbitrary linear combination of all the nj. Fhese integrals are general- 
izations of those encountered using Hylleraas basis functions, and include all that are needed to 
make energy computations on the Li atom and other four-body systems with a fully exponentially 
correlated Slater-type basis of arbitrary quantum numbers. The only quantities needed to start 
the recursion are the basic four-body integral first evaluated by Fromm and Hill, plus some eas- 
ily evaluated three-body "boundary" integrals. The computational labor in constructing integral 
sets for practical computations is less than when the integrals are generated using explicit formulas 
obtained by differentiating the basic integral with respect to its parameters. Computations are fa- 
cilitated by using a symbolic algebra program (maple) to compute array index pointers and present 
syntactically correct FORTRAN source code as output; in this way it is possible to obtain error- free 
high-speed evaluations with minimal effort. The work can be checked by verifying sum rules the 
integrals must satisfy. 

PACS numbers: 31.15.ve,31.15.vj,02.70.-c 



I. INTRODUCTION 

As long ago as 1929, Hylleraas [l[ presented 
a computation of the electronic structure of 
the He atom showing that a basis of explicitly 
correlated wavefunctions provided a far more 
efficient representation of that system than was 
available from a conventional orbital basis. What 
has come to be known as a Hylleraas atomic basis 
consists of functions, each of which is a product 
of exponentials in the electron-nuclear distances 
(often kept identical for all basis members) to 
which is appended a product of powers of both 
the electron-nuclear and the electron-electron 
distances. Hylleraas-basis computations of the 
electronic structure and properties of two-electron 
systems (i.e. the He isoelectronic series) have 
by now been successfully carried out to great 
precision by the inclusion of up to several 
thousand terms in a basis-set expansion. Repre- 
sentative results are those of Yan and Drake 0, 9 • 

An alternative to the traditional Hylleraas 
expansion is the use of basis functions that 
have correlation in the exponential, i.e. in which 
both the electron-nuclear and electron-electron 
distances appear exponentially. This type of 
basis exhibits (at modest expansion lengths) an 
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even more efficient representation of two-electron 
problems than does the original Hylleraas basis, 
and has in addition the theoretical advantage 
that, because it provides similar descriptions 
of all the particle pairs, it is also applicable 
to so-called nonadiabatic systems in which all 
three particles have comparable mass. Extensive 
computations of two-electron systems in these 
exponentially correlated bases have been reported 
by a number of investigators; representative of 
this work is a contribution by Frolov and Smith [J] . 

Calculations in the Hylleraas, the exponentially 
correlated, and other bases (e.g. containing 
logarithmic terms Q) have now been carried out, 
at least for the neutral He atom, to truly extreme 
accuracy. The situation has been summarized 
recently by Schwartz Q. 

A related line of endeavor has been to search for 
wavefunctions which yield optimum results when 
restricted to highly compact forms. Moderate 
success in this direction has been obtained using a 
basis that takes full cognizance of the asymptotic 
and other limiting behavior of the wavefunction 
0]; a greater degree of quantitative success has 
been achieved for the He isoelectronic series by 
careful optimization of four-term exponentially 
correlated functions [1, Q . 

Part of the reason for the great success with 
three-body (two-electron) problems has been that 
the necessary integrals for both Hylleraas and ex- 
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ponentially correlated functions are relatively sim- 
ple, and the organization of the integral computa- 
tions has been facilitated by the existence of re- 
cursive procedures [101 ] enabling integrals of the 
form 
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to be constructed systematically from those of 
smaller ni,n2,rii2- Here rj is the magnitude of 
and r y = |r» - rj|. 

The situation becomes drastically different for 
problems containing more than three particles. 
Historically, integrals for fully exponentially corre- 
lated wavefunctions were regarded as intractable, 
while the corresponding integrals for Hylleraas 
wavefunctions could only be evaluated by writing 
the pre-exponential powers of the inter-electron 
coordinates as spherical-harmonic expansions 
PH . Thus far, the most accurate studies of a 
four-body system, the Li atom, have used the 
Hylleraas basis set. A good survey of the current 
situation is in a review by King to which 

should be added a recent paper by Puchalski and 
Pachucki |13j ] that reports a more accurate Li 
computation than those discussed by King. 

A major advance for four- b ody systems oc- 
curred when Fromm and Hill [14j presented in 
1987 a closed formula for the basic exponentially 
correlated integral 



h{u\, U 2 ,U 3 ,Wl,W2,W 3 ) 



dxi d,Y2 dr 3 
64^3 



f ,-w 1 r 1 -w 2 r2-w 3 r 3 -u 1 r2 3 -u 2 r 13 -u 3 r 1 2 

x € - . (2) 

rir 2 r 3 r23ri3ri2 

However, the Fromm-Hill formula, though truly 
a mathematical tour de force, was difficult to 
use, even after simplifications introduced by the 
present author (l5l . Il6ll . and it was fortunate that 
in 1991 Remiddi [171 ] provided a much simpler 
formula for the basic four-body Hylleraas integral 
(that of Eq. © with the parameters itj set to 
zero) . 

While the original Fromm-Hill formula could 
be differentiated with respect to the Wi and Ui to 
introduce pre-exponential powers of the and , 
the lack of Ui dependence in the Remiddi formula 
made such an approach unavailable there. This 
difficulty was removed when Pachucki, Puchalski, 
and Remiddi [l8j published a set of recurrence 
relations enabling arbitrary increases to all the 



pre-exponential powers in four-body Hylleraas 
integrals. While Pachucki et al. indicated that the 
extension of their results to the fully exponentially 
correlated case would be "of great interest" , they 
did not consider that problem in their work. 

The present work builds upon a prelimi- 
nary study by the present author [l^] which 
presented some identities (which could be char- 
acterized as sum rules) connecting four-body 
exponentially correlated integrals with contigu- 
ous pre-exponential powers. The main result 
of the present communication is a family of 
recurrence formulas which enable construction of 
exponentially correlated integrals with arbitrary 
pre-exponential powers, starting from the basic 
integral, Eq. ([2|), and "boundary" integrals involv- 
ing fewer than four particles. It thus consitutes a 
generalization of the valuable result of Pachucki 
et al. 

While the integrals explicitly discussed in this 
paper [i.e., those represented by Eq. ©] involve 
only the interparticle distances and are therefore 
independent of the coordinates needed to describe 
the overall angular dependence of a four-body 
wavefunction, it was pointed out by Fromm and 
Hill [3] that if spherical-harmonic angular func- 
tions are included, integration over their coordi- 
nates can be carried out, leaving resultant forms 
that can be identified as cases of Eq. © . Details 
of this reduction have been addressed in previous 
work by the present author [l^jHH, so that in prin- 
ciple the technology to address P, D, ... states 
is complete. However, the evaluation of the angu- 
lar contributions to the kinetic-energy matrix ele- 
ments is complicated when expressed in terms of 
the interparticle coordinates, and there is room for 
further analysis to identify straightforward meth- 
ods for treating these states. 



II. PROBLEM FORMULATION 

The integrals that are the subject of this study 
are of the general form 
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and can be interpreted as describing the interac- 
tion of one particle (Particle 0), at the origin of 
the coordinate system, with three others (1,2,3) at 
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the respective positions r^. The integrals / there- 
fore have not only the symmetry corresponding 
to renumberings of Particles 1, 2, and 3, but also 
that resulting from rewriting Eq. |[3j to place a 
particle other than Particle at the coordinate 
origin. Specifically, the renumbering of 1-3 yields 
the identities 

fni ,ri2 ,Ti3 ,mi ,7712 ,771.3 
/ri2,ni ,n.3, m2, mi ,rri3 

(u 2 ,u 1 ,u 3 ,w 2 ,w 1 ,w 3 ) = 
J '713 ,n 2 ,ni ,7773 ,m 2 ,1111 (u 3 ,u 2 ,u 1 ,w 3 ,w 2 ,w 1 ) = 
(u 1 ,u 3 ,u 2 ,w 1 ,w 3 ,w 2 ) = 

fn2 ,ri3 ,ni ,rri2 ,m.3,mi 

(u 2 ,u 3 ,u 1 ,w 2 ,w 3 ,w 1 ) = 

/713, 711,772, 7713,7711,7712 

(u 3 ,u 1 ,u 2 ,w 3 ,w 1 ,w 2 ). (4) 

The placement of a particle other than Particle 
at the coordinate origin yields the additional rela- 
tions 

/ni,n2,ti3,mi,ra 2 ,m3 { u l, U 2 , U 3 , W\ , W 2 , W 3 ) = 
jra \ .m-2 ,773 .ni ,772 ,7713 

/m 1 ,n2,m3,ni ) m2,n 3 (^1, «2, W3, Ui, W2, U3) = 
/ni,m2,m3,mi,n2,n 3 

and the complete symmetry of the / is the 
24-element group (isomorphic with that of the 
6-j symbol) that is the direct product of the 
symmetry operations identified in Eqs. ^ and 
Notice that the parameter set (ui, u 2 ,u 3 ) 
does not have the same symmetry properties as 
(wi,w 2 ,w 3 ); the Ui relate to that form a 
triangle, while the Wi relate to that form a star. 

An important consequence of the symmetry 
relations is that it is only necessary to derive one 
key recurrence formula, which we choose to be 
that which increases the index ri\ from those of 
a reference set. Formulas for the advancement of 
all the other indices can then be obtained by an 
appeal to symmetry. 

It is convenient, following Pachucki et al, to 
define a shell of integrals as those with a common 
value of N = mi + m 2 + m 3 + m + n 2 + n 3 and 
refer to N as the shell index. We shall find that 
the key recurrence formula relates one integral in 
the shell of index 7V+1 to a number of integrals 
in shells of index N or less, so a systematic 
procedure for generating integrals in shell JV+1 
will involve the prior generation of all integrals in 
the shells with indices < N. 



Because the number of parameters and indices 
is rather large, increased compactness and clar- 
ity in the exposition can be achieved by the ju- 
dicious use of notational conventions. We there- 
fore introduce the notion of a reference index set 
ni,n 2 ,n 3 ,mi,m 2 ,m 3 and adopt the convention 
that when ambiguity will not thereby result, in- 
dices having their reference values will be omit- 
ted. We also suppress the parameters u, and Wi 
whenever possible. Thus, for example, 

(6) 

fm+l,m3 — 1 = fni,n2+l : n 3 ,m 1 ,m2,m 3 — 1 ; (7) 
/rii + l,m 3 = l = fn 1 + l,ri2,n3,m 1 ,m2,l- (8) 

The recurrence formula we shall derive is most 
directly formulated in a notation in which the 
boundary integrals entering the formula are iden- 
tified as degenerate cases of the /. Accordingly, 
using a notation introduced by Pachucki et al, we 
define 
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Note that in Eq. (J9j) , the asterisk indicates the 
presence of 47r<5(r 2 3) in place of rjg -1 exp(— 1*1^3). 
Other placements of the asterisk correspond to 
making this substitution with respect to other 
or n. The notational convention of the preceding 
paragraph also applies to these degenerate /, so, 
for example, 
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III. RECURRENCE FORMULA 

We present here the key recurrence formula, de- 
ferring its detailed derivation to Section [V] This 
formula, for / ni +i,nj,n 3 ,)ni,nu,m3i written in terms 
of the reference indices n±, n 2l n 3l mi, m 2l m 3 and 
therefore denoted simply f ni +i, takes the decep- 
tively simple form 



/rai+l — 



C\X\ + C2A2 + C3A3 
D 



(11) 
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The coefficients C±, C 2 , C 3 , and D are indepen- 
dent of the index values and are given by 

C x = ui(a&-4u|u§), (12) 

C 2 = ^2(2^3^12-^13^23), (13) 

C 3 = u 3 (2m2Mi3 - M12M23), (14) 
D = 2u 2 u 3 (u 2 1 f4 3 + U2M13 + u 3Mi2 

- M12M13M23 - ^ujujul) , (15) 

where the new quantities mj are 
I 
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The numerator quantities Xi in Eq. (fTTj) depend 
on the reference index set and on the / from shells 
of index < AT, thereby imparting the recursive 
property. The Xi have the following explicit form, 
in which the jk sum is over the two ordered pairs 
in which j and k are the members of (1,2,3) other 
than i, and 8 p is unity if p = and zero otherwise: 



X t = 2u j u k (n l + rij + 7i k + l)f - rijUk 2m f, 
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The last line of Eq. (fl"9)) contains boundary 
integrals of the type introduced at Eq. ((9]). 
As shown in Appendix B, these terms can be 
written in terms of the three-body integrals 
r ni , n2 ,„ 12 (a,/3, 7) given in Eq. Jl}. We have 

f*,n 2 ,n 3 ,m 1 ,m 2 ,m 3 — 
Tmi ,777.2+7773 — 1,712+713 — 

i(w 1 ,w 2 +w 3 ,u 2 +u 3 ), (20) 

fn 1 ,n 2 ,n 3 ,rn 1 ,m 2 ,* = 

^mi+n 2 -l,ni+m 2 -l,n 3 [w\+U 2 , Ui+W 2 , U 3 ), (21) 

and further formulas obtainable by simultaneous 
permutation of the m, Wi, mi, and m. 

The expressions given above provide a formal 
route to all / of shells with A^ > from the single 
basic N = four-body integral /oooooo = ^0 an d 
various three-body integrals T. To make this 
paper self-contained, recursive formulas for the 
r are included in Appendix B, and evaluation of 



the basic integral Iq is treated in Appendix C. 

The recursive scheme outlined above will fail 
when the quantity D or any of its permutational 
analogs are zero, a condition that occurs if any of 
the m or Wi vanish. The methods reported here 
are therefore not directly applicable to the Hyller- 
aas basis (in which all the m are zero); that case 
is more appropriately handled by the formulas of 
Pachucki et al. 



IV. NUMERICAL EVALUATION 

It is considerably more complicated than it may 
at first appear to make actual calculations based 
on the recursive process defined in the preceding 
section. Nevertheless, the recursive process turns 
out to be less cumbersome than procedures that 
depend upon the explicit evaluation of high-order 
derivatives of the basic integral presented as 
Eq. ©. 
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The applications we presently contemplate 
involve the use of basis sets that can mimic the 
ls 2 2s ground-state electronic structure of the 
Li atom, and therefore require computations at 
least as far as the shell of integrals with N = 8. 
To reach the integrals needed from the N = 8 
shell requires the evaluation of approximately 
700 integrals, and it is desirable to carry out the 
computations in a way that does not include an 
unacceptable level of organizational overhead. 

The actual approach we employed was to use 
MAPLE [22j to do the index arithmetic needed 
to write each specific instance of Eq. ([11"]) in a 
form requiring no index computations, following 
which we arranged to have these equations 
output in a form fully compliant with fortan-95 
language specifications and involving no nested 
loops. These procedures may seem to be overkill 
until it is recognized that index computations 
may require nearly an order of magnitude more 
computer time than the subsequent formation of 
the recurrence formulas. 

Grouping the FORTRAN formulas into sets 
with the same shell index, we were then able 
to carry out the recursive computations in a 
permissible order. This strategy caused the 
generation, through the N = 8 shell, of nearly 
10,000 lines of error-free code. To avoid an 
excessive accumulation of round-off error, all 
the FORTRAN computations were carried out in 
quadruple-precision floating point, and checked 
for adherence to the sum rules reported in earlier 
work [h| . The final integral values were generally 
found consistent to at least double-precision 
accuracy. We note that for the problems for 
which the methods of this paper are appropriate, 
it would not be cost-prohibitive to carry out the 
arithmetic operations with even higher-precision 
arithmetic. 



V. DERIVATION OF RECURRENCE 
FORMULA 

Following Pachucki et al. [l8|], we introduce a 
set of integrals G, of definition 



G 
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The relation between G and the integrals Y and 
/, respectively introduced at Eqs. |T]) and ([3]), is 
discussed in Appendix A; results needed here are 
Eqs. (|A9j) - l|Allj) and their permutational analogs. 

Continuing the path of Pachucki et al, we con- 
sider the following integral, which can be shown 
to vanish by application of Gauss's theorem: 
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Carrying out the operations implied by the inte- 
grand and identifying the result in terms of the G 
(a process that requires the use of identities such 
as 2ki • k 2 = k\ + k\ — kf 2 ), we reach 
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At this point it is convenient to modify Eq. ([24|) to 
a symmetry-equivalent equation by interchanging 

U2 <-> 1V2, U3 <-> W3, n 2 «-> m2, ns <-> 7713 , thereby 
obtaining 

2u 1 G2imi + M12G121111 + A ( i3Gii2iii — Gnmi 



Gqi 
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G12 



0. (25) 



We now replace the G by their equivalents in terms 
of /, using formulas from Appendix A. After mul- 
tiplying through by the factor needed to clear all 
variables from the denominators, Eq. (|25[) becomes 

2UiM 2 M3/i00000 + M12W3/010000 

+ ^13^2/001000 = Xi, (26) 

where 

Xi = 2W2U3/000000 + M3(/*ioooo — /01000*) 

+ 7i2(/*oiooo — jooioo*)- (27) 

Our next step is to apply to both sides of 
Eq. ([26]) the differentiation operator 
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n 

i=l 



_d_ 

du, 



d 
dw, 



G 



after which we define the reference index values 
to be (ni, n 2 , n 3 , J7ii, rri2, m 3 ). Looking at Eq. ([3]), 
we see that differentiation of / with respect to 
Ui (or Wi) will cause its index n, (or mj) to 
be increased by unity. Therefore, the left hand 
side of the resulting equation will contain one 
term in which the differentiations are all applied 
to the function /loooooj this term will have the 
same coefficient as the /100000 term of Eq. 
and, in terms of the reference index values, 
/100000 becomes fm+i- Similar observations 
apply to /010000 and /ooiooo- There will also be 
additional terms that result when one or more 
of the left-hand-side differentiations are applied 
to the coefficients 2u\u 2 u 3 , /112M3, or [ii 3 u 2 . We 
transpose these terms to the right hand side and 
combine them with the result of differentiating X\ . 

When V is applied to X\, we encounter differ- 
entiations of quantities such as /* 10000 ■ Keeping 
in mind that /* 10000 does not depend upon u\ but 
depends exponentially on the other Ui and W4, 
we see that X>/*ioooo will vanish unless n% = 0, 
and nonzero values of the other rii and rrii will 
result in incrementation of the non-asterisked 
indices. This lack of n\ dependence leads to the 
introduction of a factor S ni in the differentiation. 
Corresponding observations apply to the other 
terms containing asterisks. 

Based on the analysis of the preceding two para- 
graphs, the application of T> to Eq. (|26[) can be 
seen to yield the first of the three equations shown 
below. The second and third of these equations 
follow by permutation of the indices in the first 
equation. 



2uiU 2 U 3 f ni + i + /il2W3/n 2 +l + Ml3"2/n 3 + l 
Ml2W3/m+l + 2uiU 2 U 3 f n2 + i + ^23"l/ri 3 + l 
/il3W2/m+l + M23Wl/n 2 + l + 2wiW 2 "3/n 3 + l 



X 2 , 

x 3 . 

(28) 



The Xi have the values given in Eq. 1(191 



Finally, we solve the equation set, Eq. (|28j) . Ap- 
plying Cramer's Rule, we get for f ni +i- 



fni- 



1 

D 



X\ ^1 2 U 3 [l\ 3 U 2 
X 2 2uiU 2 U 3 

X 3 M23W1 2mu 2 u 3 
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with 



D = 



2mU 2 U 3 M12U3 ^13^2 
Hl 2 U 3 2uiU 2 U 3 A*23"l 
Ml3"2 M23U1 2u\U 2 U 3 



(30) 



Expanding the determinants and dividing the nu- 
merator and denominator of Eq. (j29[) by — u\, we 
obtain the expression for f ni +i shown in Eq. (|11| . 
with Ci, C 2l C 3 , and D as given in Eqs. ((12(1 -((1^ ]) . 
Wc need not exhibit solutions for /„ 2 +i or f n3 +i 
because they can be reached by permutation of the 
indices in the expression for f ril +i. 
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APPENDIX A: FOURIER 
REPRESENTATION FORMULAS 



The formulas in section [V] have forms that de- 
pend crucially on the Fourier-representation forms 
of four-body integrals of the generic type 



L = 



dridr 2 dr 3 
64^3 



x ^23(^23)^13(^13)^12(^12)^^1)^^2)^3 (7-3) 
dri dr 2 dr 3 f dqi dq 2 dq 3 dq 4 dq 5 dq 6 



64tt 3 J (2tt) 18 
X h^iq^h^q^h^iq^hi (q 4 )h 2 (q 5 )hs (q e ) 
x exp [i qi • (r 2 - r 3 ) + q 2 • (r 3 - n) 

+ q 3 • (ri - r 2 ) q 4 • Vi - q 5 • r 2 q 6 • r 3 



(Al) 

Here h(r) is a direct-space function and h T (q) is 
its Fourier transform. The transform pairs needed 
here are 



h(r) 



h(r) = 6(r), 



h T (q) 



47T 
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Performing now the integrations, which are all 
of the generic type 



J e lq r dr = (2^) 3 <5(q), 



(A4) 
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and then evaluating the integrals over q4, qs, and 
q6, we find 



L = 



dqi dq 2 dq s T T T 

h 23 {Ql) h 13 fej/liafe) 
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x hj(q23)hj(qi3)h%(qvi), 



(A5) 



where qij = |q< — qj\. 

We now insert into L as given by Eq. (| A5|) , fac- 
tors h T of the form in Eq. (|A2[) , with the result 
that L becomes equal to the integral Gmm as 
defined in Eq. (|22|) . In addition, we can insert 
the corresponding functions h into the direct-space 
form in Eq. (|A1[) . thereby also identifying L as 
/oooooo, defined in Eq. ([3]). Equating these forms 
for L, we reach 



G11 



oooooo- 



(A6) 



Next, we consider the result when we evaluate 
L taking h 2 3 and of the form in Eq. (|A3|) . 
with the other h and h T continuing as instances 
of Eq. (|A"2j) . We then have, from Eq. (|A"5|) . 
£ = Goniii/47r. Alternatively, the direct-space 
formula for this L can be identified as /*ooooo/47r, 
where the asterisk-containing / is the degenerate 
form introduced at Eq. ©. Equating these alter- 
nate forms for L, we have 



G011111 — f. 



*ooooo • 



(A7) 



Similar operations can be carried out if L is 
evaluated taking Eq. (|A3|) for /13 and h^, with 
Eq. ![A"2"]) for the other h and h T . The result is 



Giimo — /( 



00000* 
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Now, differentiating both sides of Eqs. (|A7[) and 
(jA8|) with respect to u 2 , we obtain the following 
results needed in the main text: 



G021111 — 
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/* 10000 
2u 2 ' 

/01000* 
2u 2 



(A9) 



(A10) 



Finally, we need the result of differentiating 
Eq. (|A6|) with respect to Ui: 



G 2 i 



1111 



/: 



100000 
2u± 



(All) 



Results analogous to those in Eqs. (|A9|) — (| Al 1|) can 
be obtained by simultaneous permutation of the 
first and second groups of three indices in / and 
G and the indices of u and w. 



APPENDIX B: THREE-BODY INTEGRALS 

In order to carry out the recursive process de- 
fined by Eq. (fl"9|) . we will need to evaluate integrals 
of the form introduced in Eq. (|9|). Carrying out 
the r3 integrations, two such integrals reduce to 
the three-body integrals 

_ /' dV\ dv-2 mi _i m 2 +m 3 -2 



x j,n2+"3— 2g— lUiri — (U)2+U)3)r2-(U2+U3)ri2 (Bl) 



dri dr 2 mi 

Jnl.ri2.n3. ml, m2.* — / -, ri o 7*1 



+ri2— 2 m+m2- 2 



X r 



"3 - 1 „ ~ (wi +u 2 )ri — (ill +U) 2 )r2 — «3'"12 
12 e 



(B2) 



These integrals can be respectively identified as 

-L m 1 , m 2 + 'm 3 — 1 , n 2 + '^3 — 1 (Wl, W2+W 3 ,'«2+W3), 
J- mi+n 2 — l,ni+m 2 — l,n 3 (wi+U 2 ,Ml+W2,U3), 

as shown in Eq. (|2ip of the main text. 



The asterisked / needed for the present work 
are equivalent to Y in which no more than one 
of the indices is negative (with the only negative 
value —1). The recursive methods most often 
used for evaluating T do not directly permit 
advancement of an index from —1; we also note 
that r is invariant with respect to simultaneous 
permutation of its indices and arguments. We 
may therefore identify the T needed here as falling 
into two cases: (1) r ni! „ 2ini2 with all indices 
non-negative, and (2) r_i ! „ 2j „ 12 with n 2 and ni2 
non-negative. 

For the first case, the recursive process can start 
from Toooj which by direct integration is found to 
have the value 



r 00 o(a,^,7) 



which we rewrite 



(a + /3)(a + 7 )(/? + 7 ) 



, (B3) 



r I a \ • B ooo(a,/3,7) 

r 000 (a,/3,7) = — , (B4 

a + p 

R ( a \ A)oo(o^,7) /t> pr\ 

B oo(a,P,7) = ; (B5) 



a + 7 



1 



(B6) 



(3 + 1 
We now introduce 

da) { 8(3 J { d 1 ) ' 

(B7) 
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and apply the recursive procedure of Sack, 
Roothaan, and Kolos [Toj j . leading to the follow- 
ing formulas: 



the recurrence formulas become 
1 



-l,n 2 ,n 12 



/3 + 7 



ni?i2ni2 — 71-2^12-^000 



a + P 



niT ni . 



~l~ ^2^ni ,n 2 — X,T13 ~t~ -^nm 2 tll2 



(B8) 



^niri2ni2 — ri2 ti 1 2 -^000 



1 



+ 7 



f^l-^ril — 1,712 i n 3 



. (B9) 



<W"2 + "12)! 
_)_ ^/)"l+ni2 + l 



It is a computationally stable procedure to 
evaluate first array ^4, then £?, and finally T. 



For the second case, namely the integrals 
r-i,n 2 ,ni 2 j a starting formula, again by direct in- 
tegration, is 



r n ~\ ln(q + 13)- In (a + 7) 

1 -1,0,0(0;, P, 7) = p 2 _ 7 2 ' ( BU ) 



-|- TT-ia^ — X, 712, 7X12 — 1 ^712^12 



, (B15) 



/3-7 



n 2 G 



n 2 — l,ni2 



■ Tl\ 2 G n2 ni2 — 1 ~t~ ^-Tl2tli2 



(B16) 



K n2 ri 12 — S n2 S ni2 Koo 

_ £„ 12 (1 - (5 n2 )(n 2 - 1)! 
(a + /3)™ 2 

| ^« 2 (1 - ^ 12 )(ni2 - 1)! 
(a + 7)™ 12 



(B17) 



(BIO) For j3— 7 small, it is more advisable to introduce 

x = ~(/3 + 7), 3/ = i03 - 7), to write 

ln(q + as + y) - ln(q + z - y) 
^0,0 = 7, > ( m 8J 



2y 

and to expand in powers of y. The result is 

,,2fc 



G ,o - J! 



y 



fe=0 



(2fc + l)(a + ar) 



2fc+l 



(B19) 



Differentiation of Eq. (|B19p leads to the expansion 



G„ 2 „ 12 = n 2 !ni 2 ! 



1 



ni2 — n 2 



n 2 + n 12 + 1 n 2 + n 12 + 2 



(n 2 - n 12 ) 2 + n 2 + n 12 + 2 2 

- y 



2(n 2 + n 12 +3) 
which can then be inserted into Eq. (|B15jl . 



, (B20) 



If /3 — 7 is not too small, one can proceed by a 
variant of the procedure of Sack et al. Writing 



r -i,o,o(a,/3,7) = 



Goo (a, /3, 7) 
P + l ' 

K 00 (a,f3,j) 



(B12) 



(B13) 



/3- 7 ' 
ln(a + /3) - ln(a + 7) , (B14) 



APPENDIX C: BASIC INTEGRAL I 

The integral Jo, defined in Eq. ((2j), is needed to 
start the recursive process. As discussed in [13, 
EE [l6j . the evaluation depends upon whether the 
quantity a is real, where 

2 222, 222, 222, 222 

a = Wi"W2W 3 + w 1 u 2 u 3 + w 2 UiU 3 + w s UiU 2 

...22/ 2 ■ 2 2 2 2 2\ 

+ MJjUj (w x + u 1 — w 2 — u 2 — w 3 — u 3 ) 

,22/2,2 2 2 2 2\ 

+ w 2 u 2 (w 2 + u 2 -w 1 -u 1 ~w 3 -U 3 ) 

,2 2/ 2,2 2 2 2 2\ //-ii \ 

+ w 3 u 3 (w 3 +u 3 - W 1 -Ui - W 2 -u 2 ). (CI) 
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For real a, Iq is given by 

3 



h 



1 

4^ 




where 

v(z) — sign(z) 





l-\z\ 




1 + 1*1 



7T 

12 



Li 2 



1 



2 



1 + 1*1 



(C3) 



and Li2(z) is the dilogarithm (see Formula 27.7.1 
of [HI; also Lewin [24j])- Both the logarithm 
and Li2 are multiple- valued, but, contrary to the 
orig inal formulation that required branch tracking 
[14| . Eq. (|C2p can be evaluated straightforwardly 
with all functions assigned their principal values. 



For imaginary <t, which occurs for physically rel- 
evant parameter values, Iq is obtained from 



h 



1 



2^Cl 2 ( 7 r-2tan- 1 (r i /| ( 7|) 
i=i 

Cl 2 ( 7 r-2tan- 1 ( 7 f /\a\) 



(C4) 



Here 012(6*) is the Clausen function ([23|], Formula 
27.8.1). Note that because CI2 is periodic with 
period 2ir, the presence of the arctangent does 
not cause multiple- valuedness in Eq. (|C4|) . 



The quantities Ti (i = 1,2,3) appearing in 
Eqs. (|C2| and (|C4)) are defined as follows: 



[W 2 + (Uj + U k )(Wj + Wfc)][M 2 + (Uj + W k )(Wj + ttfc)] 
Uj + Mfe + lUj + Wfc 



~ ("Uj + Ufc + Wj + Wk)(UjWj + UkWk) (C5) 



The four 7 j are 

7^ 0) = 2uiu 2 u 3 + ui/i 23 + + 

7 1 0) = ~2uiu 2 u 3 - ui^23 + + "3^12, 

72 0) = -2uiu 2 u 3 + 1*1^23 - U2M13 + "3^12, 

73 0) = -2uiu 2 u 3 + 1*1^23 + U2M13 ~ "3M12, 

(C6) 

where defined in Eqs. |[I5 )H (I5 |) . The 

jj with i =/= can be obtained from by, for 
2 = 1, the simultaneous permutation M2 <-> W2 



and U3 «-> w 3 ; for i = 2, u± <-> wi and M3 *^ W3; 

and for i = 3, u± <-> wi and u 2 w 2 . This recipe 

(i) 

produces the 7^ with a different indexing than 
in earlier work, but the value of Jo is not affected 
thereby. 



There are problems with the numerical evalu- 
ation of Iq when the parameters iWj and U{ ex- 
actly or approximately satisfy certain relation- 
ships; these situations and methods for the avoid- 
ance of numerical instability have been discussed 
elsewhere [HI, HH . 
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